package scu.maqiang.fes;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class FES1L2CPTE extends FES{
    public FES1L2CPTE(Mesh mesh, FES1L21 fes1) {
        super(mesh, 3);
        this.fes1DL21 = fes1;
    }
    /**
     * @brief 组装考虑振动与而不考虑阻尼的直接耦合法的热压电耦合整体迭代矩阵
     * @param constCoef 材料系数向量, 结构均匀各向同性，向量至少有9个分量, 分别表示
     *                  constCoef[0] - 材料密度rho
     *                  constCoef[1] - 比热容c
     *                  constCoef[2] - 热传导系数kxx
     *                  constCoef[3] - 参考温度 TRef
     *                  constCoef[4] - 热传导时间theta-迭代格式中的theta值
     *                  constCoef[5] - 弹性模量 E
     *                  constCoef[6] - 泊松比 Nu
     *                  constCoef[7] - 热模量 Beta
     *                  constCoef[8] - 热电系数 g
     *                  constCoef[9] - 压电系数 pe
     *                  constCoef[10] - 介电系数 dep
     *                  constCoef[11] - dt 迭代时间步长
     *                  constCoef[12] - alpha Newmark迭代格式参数
     * @param tpElasticStiff 组装的弹性刚度矩阵类型
     * @param tpElasticMass 组装的弹性质量矩阵类型
     * @param tpHeatStiff 组装的热刚度矩阵类型
     * @param tpHeatMass 组装的热质量矩阵类型
     * @param tpCoupled 组装的热力耦合矩阵类型
     * @param LHS 迭代左端矩阵
     * @param RHS 迭代右端矩阵
     */
    public void assembleDierctCouplingMatrix_Type2(double[] constCoef, BVPType tpElasticStiff, BVPType tpElasticMass,
                                                   BVPType tpHeatStiff, BVPType tpHeatMass, BVPType tpCoupled, SRMatrix LHS, SRMatrix RHS) {
        double rho = constCoef[0];
        double c = constCoef[1];
        double kxx = constCoef[2];
        double TRef = constCoef[3];
        double theta = constCoef[4];
        double E = constCoef[5];
        double Nu = constCoef[6];
        double Beta = constCoef[7];
        double g = constCoef[8];
        double pe = constCoef[9];
        double dep = constCoef[10];
        double dt = constCoef[11];
        double alpha = constCoef[12];

        double coefM = rho * c / TRef;
        double coefK = kxx / TRef;

        int dofHeat = mesh.getnPerElement();
        int dofElastic = mesh.getnPerElement() * mesh.getnDim();
        int dofElectric = dofHeat;
        double[][] matP = new double[1][dofHeat];
        Arrays.setAll(matP[0], i -> rho);

        double[][] ENu = new double[2][dofHeat];
        Arrays.setAll(ENu[0], i -> E);
        Arrays.setAll(ENu[1], i -> Nu);
        double[][] vecBeta = new double[1][dofHeat];
        Arrays.setAll(vecBeta[0], i -> Beta);

        double[][] matK = new double[1][dofHeat];
        Arrays.setAll(matK[0], i -> coefK);
        double[][] matM = new double[1][dofHeat];
        Arrays.setAll(matM[0], i -> coefM);
        int nv = mesh.getNv();
        double[][] HeatStiff = new double[dofHeat][dofHeat];
        double[][] HeatMass = new double[dofHeat][dofHeat];
        double[][] ElasticStiff = new double[dofElastic][dofElastic];
        double[][] ElasticMass = new double[dofElastic][dofElastic];
        double[][] HeatToElastic = new double[dofElastic][dofHeat];
        double[][] ElasticToHeat = new double[dofHeat][dofElastic];
        double[][] eleLMatrix = new double[dofHeat][dofHeat];
        double[][] eleRMatrix = new double[dofHeat][dofHeat];
        double[][] eleLMatrix6 = new double[dofElastic][dofElastic];
        double[][] eleRMatrix6 = new double[dofElastic][dofElastic];
        int[] idxHeat = new int[dofHeat];
        for(int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            int[] idxElastic = getDof(ele, nDim);
            Arrays.setAll(idxHeat, ii -> ele[ii] + nDim * nv);
            FES1L21.elementElasticStiff(ele, coordRegion, ENu, tpElasticStiff, ElasticStiff);
            FES1L21.elementMass(ele, coordRegion, matP, tpElasticMass, ElasticMass);
            //FES1L2CTE::elementDirectCouplingMatrix1(ele, coordRegion, vecBeta, tpCoupled, HeatToElastic);
            FES1L21.elementStiff(ele, coordRegion, matK, tpHeatStiff, HeatStiff);
            FES1L21.elementMass(ele, coordRegion, matM, tpHeatMass, HeatMass);

            for(int si = 0; si < dofElastic; si++) {
                for(int sj = 0; sj < dofHeat; sj++) {
                    ElasticToHeat[sj][si] = HeatToElastic[si][sj];
                    HeatToElastic[si][sj] = -HeatToElastic[si][sj];
                }
            }

            for(int si = 0; si < dofHeat; si++) {
                for(int sj = 0; sj < dofHeat; sj++) {
                    eleLMatrix[si][sj] = HeatMass[si][sj] + theta * dt * HeatStiff[si][sj];
                    eleRMatrix[si][sj] = HeatMass[si][sj] - (1 - theta) * dt * HeatStiff[si][sj];
                }
            }

            for(int si = 0; si < dofElastic; si++) {
                for(int sj = 0; sj < dofElastic; sj++) {
                    eleRMatrix6[si][sj] = 1.0 / (alpha * dt * dt) * ElasticMass[si][sj];
                    eleLMatrix6[si][sj] = ElasticStiff[si][sj] + eleRMatrix6[si][sj];
                }
            }


            LHS.assemble(idxElastic, eleLMatrix6);
            LHS.assemble(idxElastic, idxHeat, HeatToElastic);
            LHS.assemble(idxHeat, idxElastic, ElasticToHeat);
            LHS.assemble(idxHeat, eleLMatrix);

            RHS.assemble(idxElastic, eleRMatrix6);
            RHS.assemble(idxHeat, idxElastic, ElasticToHeat);
            RHS.assemble(idxHeat, eleRMatrix);
        }

    }

    public FES1L21 fes1DL21;
}
